Astronomy & Astrophysics manuscript no. turbulence paper 


© ESO 2009 


September 29, 2009 





Algorithmic comparisons of decaying, isothermal, 
supersonic turbulence 

S. Kitsionas 1 -*, C. Federrath 2 ' 3 , R. S. Klessen 1 - 2 , W. Schmidt 4 , D. J. Price 5 - 6 , L. J. Dursi 7 , M. Gritschneder 8 , 
S. Walch 8 - 9 , R. Piontek 1 , Jongsoo Kim 10 , A.-K. Jappsen 1 - 7 - 9 , P. Ciecielag 8 - 11 , and M.-M. Mac Low 12 



Astrophysikalisches Institut Potsdam, An der Sternwarte 16, D- 14482 Potsdam, Germany 

Institut fur Theoretische Astrophysik, Universitat Heidelberg, Albert-Ueberle-Str. 2, D-69120 Heidelberg, Germany 
Max-Planck-Institut fur Astronomie, Konigstuhl 17, D-69117 Heidelberg, Germany 

Lehrstuhl fur Astronomie, Institut fur Theoretische Physik und Astrophysik, Universitat Wiirzburg, Am Hubland, D-97074 
Wiirzburg, Germany 

School of Physics, University of Exeter, Stacker Road, EX4 4QL Exeter, U.K. 
CSPA, School of Mathematical Sciences, Monash University, Clayton Vic 3168, Australia 

Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON, M5S 3H8, Canada 
Universitats-Sternwarte Miinchen, Scheinerstr. 1, D-81679 Miinchen, Germany 
School of Physics & Astronomy, Cardiff University, 5 The Parade, CF24 3AA Cardiff, U.K. 

Korea Astronomy and Space Science Institute, 61-1 Hwaam-dong, Yuseong-gu, Daejeon, 305-348, Republic of Korea 
N. Copernicus Astronomical Center, Bartycka 18, 00-716 Warsaw, Poland 

American Museum of Natural History, Department of Astrophysics, Central Park West at 79th Street, New York, NY, 10024-5192, 
U.S.A. 

ABSTRACT 

Context. Simulations of astrophysical turbulence have reached such a level of sophistication that quantitative results are now starting 
to emerge. However, contradicting results have been reported in the literature with respect to the performance of the numerical 
techniques employed for its study and their relevance to the physical systems modelled. 

Aims. We aim at characterising the performance of a variety of hydrodynamics codes including different particle-based and grid-based 
techniques on the modelling of decaying supersonic turbulence. This is the first such large-scale comparison ever conducted. 
Methods. We modelled driven, compressible, supersonic, isothermal turbulence with an RMS Mach number of M lms ~4, and then let 
it decay in the absence of gravity, using runs performed with four different grid codes (ENZ0, FLASH, TVD, ZEUS) and three different 
SPH codes (GADGET, PHANTOM, VINE). We additionally analysed two calculations denoted as PHANTOM A and PHANTOM B using two 
different implementations of artificial viscosity in PHANTOM. We analysed the results of our numerical experiments using volume- 
averaged quantities like the RMS Mach number, volume- and density-weighted velocity Fourier spectrum functions, and probability 
distribution functions of density, velocity, and velocity derivatives. 

Results. Our analysis indicates that grid codes tend to be less dissipative than SPH codes, though details of the techniques used 
can make large differences in both cases. For example, the Morris & Monaghan viscosity implementation for SPH results in less 
dissipation (PHANTOM B and VINE versus GADGET and PHANTOM A). For grid codes, using a smaller diffusion parameter leads to less 
dissipation, but results in a larger bottleneck effect (our ENZO versus FLASH runs). As a general result, we find that by using a similar 
number of resolution elements N for each spatial direction means that all codes (both grid-based and particle-based) show encouraging 
similarity of all statistical quantities for isotropic supersonic turbulence on spatial scales k < N/32 (all scales resolved by more than 
32 grid cells), while scales smaller than that are significantly affected by the specific implementation of the algorithm for solving the 
equations of hydrodynamics. At comparable numerical resolution (-/V palt icies ~ W ce ii s ), the SPH runs were on average about ten times 
more computationally intensive than the grid runs, although with variations of up to a factor of ten between the different SPH runs 
and between the different grid runs. 

Conclusions. At the resolutions employed here, the ability to model supersonic to transonic flows is comparable across the various 
codes used in this study. 
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1. Introduction 

Laboratory and terrestrial flui d dynamics ar e often described 
as incompressible flow (e.g., iLesieurl Il997l) : however, astro- 
physical fluids are usually characterised b y highly compress- 
ible supersonic turbulent motion s (see e.g. Elmegre en & Scalol 
12004 IScalo & Elmegreen! |2004|) . For example, the large ob- 
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served line widths in Galactic and extragalactic molecular clouds 
and star-forming regions show direct evidence of chaotic veloc- 
ity fields with magnitudes in excess of the sound speed. This 
random motion carries enough kinetic energy to counterbalance 
and sometimes ove rcompensate for the effects of self-gravity 
in the se clouds (e.g. Ballesteros- Paredes et al.ll2~007t iBlitz et all 
2007). The intricate interplay between supersonic turbulence and 
self-gravity determines the overall dynamical evolution of these 
clouds and their observable features, such as their density struc- 
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ture , the star formation rate within them, and their lifetim es (see 
e.g. lMac Low & Klessedl2004t iMcKee & Ostrikerll2007l) . 

Despite turbulence being a universal phenomenon, it is also 
one of the least understood natural phenomena. Turbulence 
arises as a result of the nonlinear terms in the Navier-Stokes 
equat ion that governs the dyna mical behaviour of gases and flu- 
ids (lFrischlll995t lL esieu n fl997l) .A self-consistent mathematical 
formulation does not exist. Thus, analytical research mostly fo- 
cuses on finding appropriate closure equations that capture the 
bulk behaviour of the system. 

As a first approach, turbulence is characterised by two spa- 
tial scales that are connected by a self-similar cascade of kinetic 
energy that occurs over the so-called inertial range. Energy is 
injected into the system on some large scale L and dissipated 
on small scales I that are comp arable to the viscous length ^ v i sc . 
For incompressible turbulence, Kolmogorovj dl94ll) described a 
simple heuristic model based on dimensional analysis that cap- 
tures the basic behaviour of the flow surprisingly well. He as- 
sumed that turbulence driven on a large scale L forms eddies 
on that scale. These eddies interact to form slightly smaller ed- 
dies, transferring some of their energy to the smaller scale. The 
smaller eddies in turn form even smaller ones, and so on, until 
energy has cascaded all the way down to the dissipation scale 
p. ' 

L viso 

In order to maintain a steady state, equal amounts of energy 
must be transferred from each scale in the cascade to the next, 
and eventually dissipated, at a rate E - tjv^/L, where vl is the 
typica l velocity on scale L and 77 is a constant determined empir- 
ically. Kolmo g orovl d 1 94 lb assumes this rate is constant through- 
out the scales, leading to vl L 1 ^, or equivalently oc kr 1 ^ for 
wavenumbers A: oc l/L. The kinetic energy in the wavenumber 
interval [k, k + dk] is E^ n (k) oc oc L 2 ^ 3 oc k" 1 ^ and conse- 
quently the energy spectrum function E(k) - dE^ n /dk oc k~ 5 ^ 3 . 
This describes the self-similar cascade of turbulent kinetic en- 
ergy. Most of this energy resides at the top of this cascade 
near the driving scale, and the spectrum drops off steeply below 
/"vise- Because of the apparently local nature of the cascade in 
wavenumber space, the viscosity only determines the behaviour 
of the energy distribution at the bottom of the cascade below 
^visc, while the driving only determines the behaviour near the 
top of the cascade on and above L. 

Supersonic flows in highly compressible gas create strong 
density perturbations. Early att empts to understand turbu- 
lence in the interstellar medium (von Weizsackdll94"3l,ll 95 It 
IChandrasekhaHll949l) were based on insights drawn from in- 
compressible turbulence. An attempt to analytically derive the 
density spectrum a nd resulting gravitatio nal collapse criterion 
was first made by IChandrasekhd (1 1 95 1 alibi) . This work w as fol- 
lowed up by several authors, culminating in the work by Sasao 
dl973l) on density fluctuations in self-gravitating media. Larson 
(1981) qualitatively applied the basic idea of density fluctuations 
drive n by supersonic turbul ence to the problem of star forma- 
tion. (BonazzolajetaTJ (U_992|) used a re-normalization group tech- 
nique to examine how the slope of the turbulent velocity spec- 
trum could influence gravitational collapse. This approach was 
combined with low-resolution numerical models to derive an ef- 
fective adiabatic index for subsonic compressible turbulence by 
iPanis & Peraul (119981) . 

In supersonic turbulence, shock waves offer additional pos- 
sibilities for dissipation. They can transfer energy between 
widely separated scales, removing the local nature of the tur- 
bulent cascade typical of incompressible turbulence. The spec- 
trum may shift only slightly, however, as the power spectrum 
(Fourier spectrum) of a step function representative of a per- 



fect shock wave is k~ 2 . iBoldvrevl d2002l) has proposed a the- 
or y of velocity st r ucture function scaling based on the work 
of IShe & Levequd d 19941) using the assumption that dissipation 
in supersonic turbulence primarily occurs in sheet-like shocks, 
rathe r than linear filament s at the centres of vortex tubes (see 
also lBoldvrev et aU2002allbl) . Transport properties of supersonic 
turbu lent flows in the astrophysica l cont ext have been d i scusse d 
by Ide Avillez & Mac Low! J2002) and iKlessen & Liiil d2003l) 
and the fractal dimension of turbulent media by Federrath et al. 
d2009bh . 

As satisfying analytic models are rare, especially when deal- 
ing with compressible and supersonic turbulent flows, special at- 
tention is drawn to numerical approaches. A wide range of meth- 
ods are used to study turbulence, ranging fr om simulating sta- 
tistica l processes such as random walks (e.g. iMetzler & Klafterl 
I2000I). re mappin g models, or certain Hamiltonian systems 
(Isich enkol |1992|) . to hydrodynamic large-eddy simulations 
(LES). In LES only the largest spatial scales are resolved di- 
rectly using a hydrodynamic integrator. For the turbulent dynam- 
ics on smaller scales, a so-called subgrid scale (SGS) model is 
utilised. Among astrophysicists, the most often used SGS model 
is numerical dissipation, i.e. performing Implicit LES. It is not 
possible to optimise the use of closure models for astrophysical 
turbulence through comparisons with laboratory experiments. 
Therefore, the representation of the SGS behaviour provided by 
numerical dissipation must be s ufficient, and inde ed provides a 
reasonably good approximation (Benzi et al. 20081). 

In the current study we focus on comparing different Implicit 
LES schemes. Our goal is to assess the applicability of differ- 
ent numerical schemes to the modelling of supersonic turbu- 
lent flows, and to compare their validity and accuracy in the 
astrophysical context. To keep this comparison simple, we fo- 
cus on purely hydrodynamic turbulence in isothermal gaseous 
media in regions with periodic boundaries, and study the decay 
of fully developed turbulence. We follow the dissipation of ki- 
netic energy due to the numerical viscosity intrinsic to any nu- 
merical scheme, and characterise the turbulent velocity field us- 
ing volume- and density-weighted velocity power spectral and 
probability distribution functions of density, velocity and veloc- 
ity derivatives. We remind the reader that in supersonic turbu- 
lence energy is not only dissipated below by the action of ar- 
tificial viscosity on the smallest scale eddies, but also in shocks. 
In most of the codes employed here artificial viscosity is nec- 
essary also for the modelling of shocks. The main aim of our 
comparisons is to characterise the role of artificial viscosity in 
dissipating energy below ff? rather than the use of artificial vis- 
cosity in the modelling of shocks. A discussion on the shock 
capturing ability of the codes used in this study is provided in 
§|5](for a comprehensive such comparison, see also lTasker et al.l 
2008). One of the fundamental questions we want to address is 
at which numerical resolution are different numerical schemes 
capable of modelling supersonic turbulence adequately. 

This is the first such comparative study; there has been no 
coherent comparison of the various hydrodynamic codes used in 
the literature for the study of supersonic turbulenc^- In spite 
of the fact that results from different codes appear to contra- 
dict each other and lead to different interpretations of the role of 



1 Volume-weighted velocity power spectra are often re fered to as ki - 
netic energy spectra for incompressible turbulence (e.g., Frisch 1995). 
However, for compressible turbulence the kinetic energy is proportional 
to the density-weighted velocity power spectrum. 

2 A limited comparis on of two codes was presented for self - 
gravitating turbulence in Kl essen et lEI fcOOCh and lHeitsch etafl fcOOll) . 
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tu rbulence in astrophysics - e.g. th e hy drodynamic simulatio ns 
of iBallesteros-Paredes etail d2006l) and IPadoan et af] d2007l) in 
which different power-law slopes are obtained from the velocity 
power spectra leading to different interpretations of the role of 
turbulence on cloud fragmentation and the resulting core mass 
function - it has not been properly checked whether at least 
some of these differences are due to differences in the numer- 
ical schemes employed. We perform and analyse here a first set 
of low-resolution calculations aiming at extending our investi- 
gations in the future to higher resolution simulations, achieved 
either directly, or by using a daptive resolution techniques like 
Adap tive Mesh Refinement (Krit suk et al.ll2006t ISchmidt et al.l 
2009), or Particle Splitting dKitsionas & Whitworthl]2002l) . or 
by using Subgrid Scale Models (ISchmidt et al. 2006b c), and/or 
combinations of the above. It should be emphasised that the typ- 
ical number of resolution elements used here for the SPH cal- 
culations is quite large (number of particles N = 215 3 ) com- 
pared to the typical number of particles used in existing studies 
of supers onic turbulence and cloud fra gmentation in the litera- 
ture (e.g. Ball esteros-Paredes et al. 120061) . On the other hand, the 
number of resolution elements used for the grid codes presented 
here is rather small (number of grid cells N = 256 3 ) compared 
to what is the current state - of-the-art resolution for such stud- 
ies (e.g. iRritsuk et al.ll2007t IPadoan et al.ll2007t IFederrath et all 
I2009at iLemaster & Stonell2009t ISchmidt et al.ll2009l) .~ 

The structure of our current study is as follows: in § [2] we 
describe the setup of the experiments conducted as well as list 
the most important features of the codes used. In § [3] we re- 
view the statistical measures used in this paper for the analysis 
of supersonic turbulence. In §|4j we present the initial conditions 
employed for the decay simulations, while in § [5]we discuss the 
results of the decay experiments and the comparison of the per- 
formance of the various codes. We summarise the computational 
efficiency of the codes and runs in § [6] and present our conclu- 
sions in § |7] 

2. Experimental setup and numerical codes 

Our aim is to study the decay of supersonic hydrodynamic tur- 
bulence using different grid- and particle-based code£j and com- 
pare the performance of the codes in this experiment. Therefore, 
we need a turbulent gas distribution that will serve as an initial 
condition for all codes. For simplicity, the turbulent initial condi- 
tions are produced with one of the particle/SPH codes, and then 
the particle code data is interpolated onto a grid. This grid, in 
turn, provides the initial conditions for the grid-based codes. 

The turbulent g as distribution is created with GADGET 
dSpringel et al.ll200 lb . We start with a box of side L = 0.29 pc. 
Inside this box, 21 5 3 = 9, 938, 375 particles were distributed 
homogeneously representing a static, uniform, isothermal gas 
with temperature T = 1 1 .4 K (corresponding to a sound speed 
c s = 0.2kms _1 ), and mass M = 120 M . We impose a turbu- 
lent veloc i ty fiel d within our box using the driving scheme of 
Mac Lowl d 1999b. as it has been implemented for GADGET by 
Jappsen et al.l d2005l) . Turbulence is driven on large scales (at 
wavenumbers between k — 1 and k = 2), aiming at an RMS 



3 We use the term particle code as the generic antonym of grid 
code. In general, a particle code is a numerical scheme that uses 
sampling points that are not fixed in space but rather move, re- 
sembling in this respect the properties of fluid particles. In particu- 
lar, all particle codes used here are different implementations of the 
Smoothed Particle Hydrodyn amic s (SPH) tech nique first introduced by 
iGingold & Monaghanl i 19771) and lLucvl d!977t) . 



Mach number equal to M lms ~ 3.9. The RMS velocity has then 
reached V rms = M rms c s ~ 0.78 km s~'. 

We drive turbulence for a few code unit times. The code unit 
for time is arbitrarily chosen to be equal to the free-fall time 
fff = ^3n/(32Gpo), of the initial homogeneous and static gas 

distribution. Using po = M/L 3 , we obtain fff ~ 0.1 15 Myr. This 
does not imply that the gas in our box is or will become self- 
gravitating. A useful time unit for the comparison of non-self- 
gravitating turbulent flows is the turbulent crossing time ? cross = 
L/(2V rms ) = L/(2M rms c s ) ~ 0.1 82 Myr, which is defined as the 
time it takes for a typical turbulent fluctuation to cross half of 
the computational domain (d efined equivalently to Kritsu k~etail 
120071: IFederrath et alJl2009M> . 

In § 14.11 we discuss for how long the turbulence needs to be 
driven to reach a statistical steady state. After driving has fin- 
ished, the density and the velocity of the GADGET-particle distri- 
butions are interpolated onto a grid with 256 3 cells. We use this 
grid as the initial condition for following the decay with the grid 
codes. The particle distribution at the end of the driving phase 
provides the initial conditions to the SPH codes. 

We follow the decay of turbulence for about six turbulent 
crossing times f cross . The self-gravity of the gas is kept switched 
off at all times and the influence of magnetic fields is ne- 
glected, i.e. we follow the hydrodynamic decay only. An isother- 
mal equation of state is used. Periodic boundary conditions are 
adopted. 

Snapshots were taken at 0.0, 0.06, 0.31, 0.62, 3.1, and 
6.2 f cross . For the SPH codes, the particle distribution is inter- 
polated onto a grid at each snapshot. The grid codes naturally 
provide their density and velocity fields on grids. From these 
grids, we then calculate spatially averaged quantities like the 
RMS Mach number, velocity power spectra, and probability dis- 
tribution functions of several quantities including the density, the 
velocity, its derivatives, and combinations of density and veloc- 
ity. 

In this paper, we include turbulence decay experiments using 
the SPH codes GADGET, PHANTOM (runs: A, B), and VINE, as well 
as the grid codes ENZO, FLASH, TVD, and ZEUS. In the following 
paragraphs, the general features of the codes used in this work 
are listed briefly. For more details on each of the codes, please 
refer to the references given below. All codes were run in paral- 
lel. The parallelisation architectures that the codes were run on, 
the total number of CPUs used, and the number of CPU hours 
consumed are listed in §|7]for each of the runs studied here. 

Codes that are used to perform simulations of supersonic tur- 
bulence are chosen for their performance in the highly compress- 
ible regime typical of astrophysical flows. Thus, codes with the 
ability to capture accurately the sharp discontinuities and high 
density ratios that result in supersonic turbulence are preferred 
for this study. 

Grid-based codes for supersonic flows are often based on fi- 
nite ^yohimeRiemann-solvers using the Godunov scheme (see 
e.g. lTorol fl997). By their conservative nature, these codes main- 
tain correct shock speeds, and, because of the Riemann-solver 
approach in calculating fluxes, they maintain very sharp discon- 
tinuities across a shock (typically within a few zones). These 
methods are often implemented in a dimensionally split way 
dStrang|ll968l) . 

A very common such method used for astrophysical 
flows is the Piecewis e Para bolic Method (PPM), described in 
IColella & Woodward! ([1984), which is formally third-order ac- 
curate in space for smooth flows but switches to linear or- 
der to maintain step-function-like shocks or contact discontinu- 
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ities near such features. Although not typically crucial, a small 
amount of artificial viscosity is impleme nted to ensure no or little 
oscillations behind shocks. Bot h ENZO ( Norman & Bryan 1999; 
O'Shea etal1l2004l) and FLASH dFrvxell et al.ll2000t iDubev et al l 
2008), used in this paper, are of this type. ENZO was used here 



in version 1.0.1, and FLASH was used in version 2.5. FLASH in 
particu lar has been exten sively tested again st laboratory experi- 
ments (|Calder et aU2002l) and other codes (iDimonte et al.l2004t 
iHeitmann et al.ll2005b .~For the FL ASH run of this work, a PPM 
diffusion parameter of K = 0.1 (IColella & Woodward! [19841) 
was used, whereas for the ENZO run the PPM diffusion param- 
eter was set to K = 0.0. This provides an additional test of the 
effects of artificial shock diffusion on the results of grid codes 
using the PPM. The effects of varying the PPM diffusion pa- 
rameter has been investigated before: setting K = 0.0 provides 
less dissipation, but produces stronger po st-shock oscillations 
and a more pronounc ed bottleneck-effect dKritsuk et al.l l2007t 
iFederrath et"aill2009al) . 

Other methods that attempt to maintain flatness of the solu- 
tion behind shocks include Total Variation Diminishing (TVD) 
methods, which impose restrictions on the reconstruction of flow 
variables to ensure that the total variation of variables strictly 
diminishes over time. The TVD code used here employs an ap- 
proximate (second-order-accurate), Roe-type (upwind) Riemann 
solver and TVD interpolation to maintain sharp shocks and 
smooth flow behind the shocks. Recipes for building the isother- 
mal MHD code based on the TVD scheme are presented in 
iKim et all d 19991) . For the turbulence comparisons of this project, 
we have used an isoth ermal hydrodynamic vers ion of the code. 

The ZEUS code (IStone & Norman 1992a,b) is a second- 
order accurate code using the Ivan Leerl ( 1 19771) monotonic ad- 
vection algorithm. It resolves shocks using artificial viscosity, 
and does not include explicit techniques for keeping shocks 
sharp. A staggered mesh approach is adopted: the velocity is 
stored at cell interfaces, and density and energy at cell cen- 
tres. This version of ZEUS is however different from the offi- 
cial version of ZEUS -MP, wh ich was employ e d for instance by 
IVernaleo & Revnoldsl d200fj|) and lHaves etal] d2006l) , in that the 
version used here was parallelised by R. Piontek. 

An entirely different method is employed in Smoothed 
Particle Hydrodynamics (SPH) codes (e.g. iLucvl 119771: 



iGingold & Monaghadl977tlMonaghadl98all992ll20()5T) . SPH 
codes do not involve a grid at all, but track (fixed)-mass pack- 
ets of fluid in a Lagrangian sense. While this approach to fluid 
dynamics requires t he use of an artificial viscosi ty to maintain 
shock structure (e.g. lMonaghan & Ging old 1983), it has the ad- 
vantage for highly compressible flows that resolution automati- 
cally increases in high-density regions, as the Lagrangian fluid 
packets follow the mass flow. Three SPH codes, namely GADGET, 
PHANTOM (runs: A, B), and VINE, are used in this work. More de- 
tails of these codes are given below. 

GADG ET is an MPI parallel tree-SPH and N-body code de- 
signed by Spri ngel et al.l d2001l) . We here used GADGET version 
1. The code uses individual and adaptive timesteps for all par- 
ticles, and it combines this with a scheme for dynamic tree 
update. As a time integrator, it uses a variant of the leapfrog 
integrator, involving an explicit predictor step. The smoothing 
lengt hs are derived according to t he commonly used M 4 ker- 
nel dMonaghan & Lattonziolll985t ISpringel et all 1200 ll) . In the 
GADGET run performed here, we set the number of neighbours 
of each SPH particle to 40 + 5 neighbours. The influence of 
changing the number of SPH neighbours A^ ne igh has been inves- 
tigated bv lAttwood et al.l (120071) . They argued that using a fixed 
number of neighbours prohibiting any variation in A^ ne i e h may re- 



duce numerical dissipation. Commerco n et al.l ((2008) also stud- 
ied the influence of changing the number of SPH neighbours 
in simulations of gravitational fragmentation. They find a weak 
dependency of their results on the number of neighbours indi- 
cating that increasing iVneigh speeds up gravitational fragmen- 
tation slightly. All SPH codes used here employed roughly the 
same number of SPH neighbours (see below). A number of about 
Mieigh ~50 ± 5 SPH neighbours is the typical setup for most SPH 
calculations reported in the literature. Thus, our comparison of 
SPH runs with each other should not be systematically affected 
by our choice of SPH neighbours. However, varying the num- 
ber of neighbours should be investigated in a detailed systematic 
study of supersonic turbulence in the future. 

PHANTOM is a low-memory, highly efficient SPH code writ- 
ten especially for studying non-self-gravitating problems. The 
code is made very efficient by using a simple neighbour find- 
ing scheme based on a fixed grid and linked lists of particles. In 
particular, it uses an rj = 1 .2 term in calculating the smoothing 
length h through h = rj(m/p)' 3 . This corresponds to about 58 
SPH neighbours in a uniform density distribution, though it is the 
h — p relation (to a tolerance of 10~ 4 ) that provides the deciding 
criterion, not the neighbour number. The c ode implements the 
full "g rad- h" SPH formulation develop ed bv lPrice & M onaghan 
d2004b and iPrice & Monaghanl (120071) . whereby the smoothing 
length and density are mutually dependent and iterated self- 
consistently, resulting in exact conservation of momentum, en- 
ergy an d entropy in the S PH equations. Shocks are treated us- 
ing the Monaghan ( 199 7D f ormul ation of artificial viscosity in 
SPH as described in iPricd d2004l) . though modified slightly in 
PHANTOM to allow a more efficient calculation. Timestepping is 
performed using a Kick-Drift- Kick leapfrog integrator. The stan- 
dard PHANTOM run is labeled as PHANTOM A. We conduct a sec- 
ond run with PHANTOM, labeled PHANTOM B, in which dissipation 
is reduced away from shocks us ing the viscosity switch proposed 
bv lMorris & Monaghanl d 1997b . 

VINE dWetzstein et al] 120081: iNelson et all [2008; 
iGritschneder et all 12009b is an OpenMP parallel, tree-SPH 
and A^-body code. The code scales linearly up to a high num- 
ber of processors, and is designed for a combined usage of 
generic CPUs and the special purpose hardware GRAPE. Time 
integration is performed here with a Drift-Kick-Drift leapfrog 
integrator, which allows for individual particle timesteps. The 
smoothing lengths are deriv ed according to the M4 kernel 
dMonaghan & Lattanziolll985l) . In the VINE simulations, we set 
the number of neighbours of each SPH particle to 50 + 5. Shocks 
are treated with the time dependent artifi c ial vis cosity prescrip- 
tion introduced by Morris & Monaghan] d 1997b (i.e. similar to 
the PHANTOM B run). 



3. Methods for the analysis of statistical measures 

3.1. SPH interpolation onto a grid 

For interpolating the density and the velocity distribution of the 
SPH particle s onto a grid, we use the generic SPH interpolation 
formula (e. g. iMonaghanll 1 992b 



V Pi h] \ 



■Ti\ 



hi 



(1) 



where p, is the density, m, is the mass, and h\ is the smoothing 
length of the z'-th particle. The vector r is the position vector 
to the centre of each grid cell, A, is the particle quantity to be 
interpolated to the grid (density and velocity for our purposes), 
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and W(s) is the kernel function used for smoothing the particle 
mass in space to derive an SPH density. The M4 kernel, which 
is based on spline functions (Monaghan & Lattanzio 1985), 



l I 1 - 3s 2 /2 + 3s 3 /4, 
W(s) = - I (2 - s) 3 /4, 
71 0, 



< s< 1, 

1 < s < 2, 
s>2, 



3.3. Density-weighted velocity power spectrum 

For density-weighted velocity power spectra, we substitute 

V mw = (P/P0)' /2 V , 



(8) 



(2) 



is used for the interpolation. For each grid point, the above sum- 
mation is over a limited number of neighbouring SPH particles 
due to the compact support of the kernel function, which van- 
ishes for |r - r,| > 2h;. 



3.2. Volume-weighted velocity power spectrum 

The velocity power spectrum is calculated as follows. For each 
velocity component we take the Fourier transform of the velocity 
fieldv = (vi, V2, V3). We denote these Fourier transforms as v' = 
(v'j ,v' 2 ,v' 3 ). Using these definitions the volume- weighted velocity 
power spectrum is defined as 



where po is the mean density in the cube. Then, the above pro- 
cess for the calculation of velocity power spectra is repeated. 
The density- weighted velocity power spectrum is defined as 



(9) 



E(k) = \ (v' • v') , 



(3) 



where v' = (y' v v' 2 ,v'^) is the complex conjugate of the trans- 
formed velocity. A wavenumber mapping, k = (£4, £2, £3), is 
applied on the cells of the E(k)-c\xh&, with each ki, ki, and ki 
ranging from -N/2 to (N/2) - 1. The compressible (longitudi- 
nal) part of the velocity power spectrum is calculated as 



in analogy to equation ||3}. The density-weighted dissipation 
rate is computed in analogy to equation ©, from the density- 
weighted velocity power spectrum E mw (k). 

3.4. Probability distribution functions 

Using all the cells in our grid, we obtain the cumulative distri- 
butions, F, of the following quantities: logarithm of the density, 
the three velocity components v, with i = 1, 2, 3, the logarithm 
of the trace free rate of strain, |5"|, and the logarithm of vortic- 
ity, u> = |V x v|. To obtain these distributions, the corresponding 
quantities are binned linearly. From the cumulative distribution, 
F, we derive the probability distribution function (PDF) of quan- 
tity A at each bin n by computing 



PDF„(A) 



F„(A)-F n - l {A) 



(10) 



E com (k) = E(k) 



(v' • k) (v' • k) 



(k ■ k) (v' ■ v') 
and the solenoidal (transverse) part as 
E soi (k) = E(k) - E com (k) . 



(4) 



(5) 



For the trace free rate of strain, the spatial derivatives of the 
velocity components are first calculated as 



'<•./ 



dXj 



(11) 



The (symmetric) strain tensor has components, S ij = 0.5(vy + 
vjj). The rate of strain is then 



For each wavenumber k = |k|, we collect from the E(k)-cube 2 V 1 V 
all cells lying at distances in the [k, k+l] intervafl The mean I" I ~ ^ 2-i 2-i V 'i 
of the £(A:)-values of these cells is normalised by the area of ' J 

the sphere element with radius k + 0.5. This gives the volume- 
weighted velocity power spectrum E(k). The above process is 
repeated for E so \(k) and E com (k). 

To calculate the volume-weighted dissipation rate, we esti- 
mate at each snapshot the integral over the volume-weighted ve- 
locity power spectrum. This is formally given as 



(12) 



The trace of the strain tensor is d = -S u, so that the trace free 
strain tensor has components S*j = Sij - dijd/3, where <?y is 
the Rronecker unit function. The trace free rate of strain then 
becomes 



(13) 



1 ■> 
- V 2 

^ rms 



E(k)dk . 



(6) 



and is estimated numerically. The dissipation rate is then given 
as the rate of change of this integral as the decay proceeds. 

The sonic scale, / s = 2nk s , separates supersonic turbulent 
flows on large scales (small k) from subsonic turbulent flows on 
smaller scales (high k). We estimate the sonic scale by solving 
the following equation, 



E(k) dk . 



(7) 



2nlk 



implicitly for the sonic wavenumber k s 



4 Distances are measured with respect to the (0, 0, 0) cell, i.e. the cen- 
tral cell. 



\S*\ gives the rate at which a fluid element is deformed without 
changing its volume, e.g. by the act of a shear flow. 



4. Driving phase 

4. 1 . Driving time 

We simulat e driven turbulence u sing GADGET following the pre- 
scription in iJappsen et alj d2005l) . The driving phase starts with 
an initially homogeneous density distribution with po = 3.3 x 
KT 19 gcirr 3 at rest. As turbulence gets driven, the RMS Mach 
number increases with time. It levels off at a value of M lms ~ 3.9 
after about 4 % (see Fig. [TJ. This time corresponds to roughly 
2.5 crossing times, / C10SS . We start the decay experiments at this 
time, i.e. after turbulence has been driven for 2.5 f C ross- The par- 
ticle distribution obtained at this time is interpolated onto a grid. 
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This grid is used as the initial condition for the decay simula- 
tions using the grid codes. For the decay simulations using the 
SPH codes, we directly use the particle distribution obtained 
with GADGET at the same time (f = 2.5 fcross)- F° r the follow- 
ing calculations the time was reset to zero, and the turbulence 
decay was followed over roughly six crossing times with each of 
the grid and SPH codes. 

Turbulence has been established when the decay experi- 
ments start. This is shown in Figure [2] where we plot the veloc- 
ity power spectra, calculated as explained in § |3.2| and §[33] for 
the volume- and density-weighted spectra, respectively, of four 
different snapshots taken along the driving phase. On the left 
panels, volume-weighted spectra are plotted whereas density- 
weighted spectra are shown on the right panels. Driving for 
1.2 fcross (black lines) was not sufficient to produce a statistically 
fully established turbulent flow, as there was not enough time for 
turbulence (driven on large scales) to cascade down to the small- 
est spatial scales. However, the turbulence is fully established af- 
ter about 2.5 fcross: there is no significant variation of the velocity 
power spectra when we attempted to drive turbulence for longer 
times (cf. the red, green and blue lines corresponding to driv- 
ing times of 2.5, 3.1, and 3.7f cr0 ss, respectively). We conclude 
that starting the decay with the gas distribution obtained after 
2.5 fcross of dr i ving i s a reasonable choice of initi al conditions- 
Sch midt et all (120091) . IFederrath et al.1 d2009bl) . and lGlover etal] 
d2009t) also conclude that after about 2 f cr0 ss, supersonic turbu- 
lence has established a statistical invariant state. However, sig- 
nificant statistical flu ctuations from snapshot to snapshot remain 
dFederrath et al.ll2b09al) . which explains the slight changes visi- 
ble in the velocity power spectra at f = 2.5, 3.1, and 3.7 f cr0 ss in 
Figure [2] The variations seen in Figure [2] for f > 2.5 f cr0 ss are at 
most in the order of the typical snapshot-to-snapshot variations 
intro duced by intermitte nt fluctuations, i.e., less than 10% (see 
also lKritsuk et al.ll2.QQ7l Fig. 1). The initial conditions used for 
our code comparison therefore constitute a statistically fully es- 
tablished supersonic turbulent density and velocity distribution. 



4.2. The result of driving: initial conditions for the decay 
experiments 

In this section, we present the velocity power spectra of the ini- 
tial conditions used for the decay experiments. These initial con- 
ditions have been produced with GADGET using the turbulence 
driving routine deve loped by Mac Lowl (1 19991) . and employed in 
iJappsen et al.l (120051) . They present the state of the system after 
2.5 fcross of driving (see § 14. j} . On the left panel of Fig. [3] we 
plot the volume-weighted velocity spectrum, with the density- 
weighted velocity spectrum shown on the right panel. The spec- 
tra were compensated with power-law slopes of 2.20 (left panel - 
volume-weighted case) and 1 .67 (right panel - density-weighted 
case). 



4.2.1 . Volume-weighted velocity power spectrum of the initial 
conditions 

From the volume-weighted velocity power spectrum computed 
with GADGET we derive a slope of about 2.2, which is obtained 
in the wavenumber range A<k< 12. If any inertial -range scaling 
could be inferred at all due to our limited numerical resolution, 
it may only exist in the close vicinity of k ~ 8. 

In the presence of a bottleneck e ffect (e.g.. | Dobler et all 

20031: iHaugen & BrandenburgJ l2004t ISchmidt et al.l l2006al) . 
Porter & Woodward! (1 19941) argued that the bottleneck affects 



all scales up t o k ~ N/32 = 8 for grids of size N = 256. 
ISchmidt et al.l (l2006al) suggested that the bottleneck peaks at 
k ~ N/10 ~ 26. These authors also argued that, in codes 
showing no bottleneck, numerical dissipation will start acting 
at wavenumbers smaller than k ~ N/10. Since our initial con- 
ditions do not seem to exhibit a bottleneck, dissipation will cer- 
tainly start at k < 26. Since a power law is established for scales 
4 < k < 12, and since this power law breaks down at k > 12, we ar- 
gue that dissipation did not play a significant role for wavenum- 
ber s k < 12 in the spectr a of th e initial condi ti ons. 

iPadoan et all d2007l) and iKritsuk et all d2007l) performed 
high-resolution hydrodynamic simulations of driven turbulence 
using ENZO at resolutions of 1024 3 grid cells. They obtained a 
significantly sh allower slope of 1 .9 for t he volume- weighted ve- 
locity spectrum. Federrath et al. (2009a) showed that driving tur- 
bulence with the FLASH code at resolutions of 51 2 3 and 1024 3 
grid ce lls r esults in a slope con sistent with those of Padoai Tet al.l 
d2007l) and IKritsuk et al.l d2007l) . In contrast, at the resolution of 
256 3 ( as used here), a steeper slope of the order of 2.1 was ob- 
tained dFederrath et al.ll2b09al) . whic h is consistent with our re- 
sult for the volume-weighted spectra. IFederrath et alj d2009al) di- 
rectly demonstrated the steepening of the velocity spectrum with 
decreasing numerical resolution (cf. their Fig. C.l). However, 
they also find that the slope of the velocity spectrum is converged 
to within less than 3% going from 51 2 3 to 1024 3 in their numer- 
ical experiments. Considering the low resolution in the present 
study, it is therefore not surprising that we find a slope of about 
2.2 for the volume-weighted velocity power spectrum. This re- 
sult however can also be taken as an indication that the initial 
conditions used for our comparison experiments did not strongly 
depend on the method by which these initial conditions have 
been produced (i.e. GADGET). 



4.2.2. Density-weighted velocity power spectrum of the initial 
conditions 

From the density-w eighted velocity pow er spectrum we obtain a 
scaling close to the iKolmogorovl dl94ll) scaling with a slope of 
about 5/3. This slope is found within 5 < k < 20. We have used 



(p/on) 1 /2 v(see §l3.3l l instea d ofv m w = (p/po) 1/3 v, which 



1/3, 



was used by IKritsuk etal] (120071). iKowal & Lazarianl d2007), 
ISchmidt et atld2008l) . and lFederrath et al.l(l2009al) . The(p/pn) 1/2 
weights correspond to a quantity that has physical reference to 
kinetic energy (l/2)p|v| 2 , while the (p/po) 1 ^ 3 weights corre- 
spond to a co nstant kinetic energ y dissipation rate within the 
inertial range (Krit suk et al.1 12007). For GADGET only, we addi- 
tionally compute energy spectra with the (p/po) 1 ^ 3 weights. The 
comparison of the (p/po) 1 ^ 2 to the (p/po) 1 '' 3 weights is shown in 
Figure |4] (left panel). We find a ste eper slope of about 1.8 for 
the (p/po) 1/3 weights. The fact that IKritsuk et all (|2007|) obtain 
Kolmogorov-type scaling using the (p/po) 1 ^ 3 weights is a conse- 
quence of their volume-weighted velocity power spectrum being 
shallower than ours with slopes of 1.9 and 2.2, respectively. We 
find this steeper slope of about 2.2 due to our limited numer- 
ical resolution as discussed in the previous section. Therefore, 
the fact that our density-weighted spectra show scaling close to 
Kolmogorov scaling is a result of the rather small numerical res- 
olution adopted in this comparison. Clearly, the slopes of the 
density-weighted spectra depend not only on the velocity statis- 
tics, but also on the convolution of density and velocity statistics. 

However, we argue that the density information should be 
taken into consideration in the statistical analysis of compress- 
ible turbulence, as most of the mass ends up in small volumes 
through shocks. This fact is neglected by statistical measures that 
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take into account volume-weighted velocities only (for instance, 
some models of the mass distribution of molecular cloud cores 
and stars a re based on the volume- we ighted velocity power spec- 
trum, e.g., Pad oan & Nordlundl2002l) . 

4.2.3. The effective SPH resolution 

We would like to comment here on the Lagrangian nature of 
the SPH method. In Figure [3] there is a rise in power on scales 
k > 100, particularly prominent in the density-weighted case 
(right panel). This is a consequence of the adaptivity in reso- 
lution that is intrinsic to the SPH method: as the SPH particles 
move with the flow, the build-up of high densities is accompa- 
nied by an increase in the number of particles (sampling-points) 
for a fixed volume element. Therefore, as high densities build up 
on small scales due to the shocks developed in supersonic turbu- 
lence, the SPH particle concentration increases on these scales. 
In particular, in driven turbulence the effective SPH resolution 
in high-density regions eventually becomes superior to the reso- 
lution initially employed. Hence, the extra power developed on 
scales k > 100 is a result of the interpolation of the SPH parti- 
cle distribution onto a grid of resolution lower than the effective 
SPH resolutioifl When the SPH particle distribution is interpo- 
lated onto a larger grid (N = 5 12 3 ) this rise in power is no longer 
observed (right panel of Fig. @J. In other words, because of the 
finite extent of the grid we used for our comparison experiments 
(N = 256 3 ), all SPH information on scales smaller than k max gets 
interpolated into the k max bin. This appears as a rise in the power 
on the smallest scales of our grid. 

5. Results of the turbulence decay code 
comparison 

5.1. Volume-weighted and density-weighted velocity power 
spectra 

Figures [5] and [6] present the time evolution of the volume- and 
density-weighted velocity power spectra obtained with the var- 
ious codes employed in this work. Data from the following 
snapshots are shown: initial conditions at t — 0.0f C ross (black 
lines), t = 0.06f closs (red lines), t = 0.31f cros s (green lines), 
t = 0.62f cross (blue lines), t = 3.1 /cross (cyan lines), and 
t = 6.2 fcross (magenta lines). Figures [7] and [8] show the volume- 
and density-weighted velocity power spectra of each code plot- 
ted on top of each other in a single plot as a function of time, 
so that the spectra obtained for each snapshot can be directly 
compared across all codes for each time. 

The grid codes dissipate the power produced on small spa- 
tial scales (k > 100) of the initial conditions faster than the 
SPH codes. This is a result of the SPH interpolation onto the 
grid (see § 14.2.31) . Due to their Lagrangian nature the SPH 
runs (GADGET, PHANTOM A, PHANTOM B, and VINE), maintained 
power at k > 20 for a longer time. All grid codes have lost 
slightly more power for k > 20 than the SPH codes, immedi- 
ately after the decay simulations start (i.e. from the first snapshot 
at 0.06 fcross)- The differences seen at early times on the small 
scales of the SPH power spectra is a result of the different meth- 
ods that each of the SPH codes adopt for calculating particle 
smoothing lengths and/or the use of different smoothing kernels, 
and different implementations of artificial viscosity. 

3 Note that the total number of sampling points for the SPH runs 
(215 3 particles) was smaller than the number of sampling points em- 
ployed for our grids (256 3 grid cells). 



The volume-weighted velocity power spectra obtained with 
the SPH codes and with the grid code ZEUS are quite similar 
for k < 20 at t — 0.06 f cro ss, with power law slopes of about 
2.2. ENZ0, FLASH, and TVD have a slightly shallower slope of 
about 2.1. This slope a grees with the low-resolution models in 
iFederrafh etal] (l2009al) using 256 3 grid cells. At t = 0.31, and 
0.62 f cross , also ZEUS develops a slightly shallower slope that 
roughly agrees with the slopes obtained using the other grid 
codes. The slopes have droped to about 1.95 and 1.9 at t — 0.31, 
and 0.62 f cr0 ss, respectively. The wavenumber ranges over which 
these slopes are maintained are slightly smaller for ZEUS than for 
TVD and FLASH (up to k ~ 12), while ENZO maintains the slopes 
up to k ~ 18. The SPH codes again have a slightly steeper slope 
by about 0. 1 than the grid codes, and the wavenumber range over 
which this slope is maintained is comparable with the range ob- 
tained using the ZEUS code. 

The density-weighted velocity power spectra are shallower 
than the volume-weighted spectra for all codes with slopes 
of about 1.6 at t = 0.06f crO ss- This much shallower slope is 
a result of the low resolution of our numerical experiments 
as discussed in § 14.2.21 Similar to the results obtained from 
the volume-weighted spectra, all grid codes dissipate the ini- 
tial power on scales k > 20 faster than the SPH codes with 
ZEUS having dissipated most. However, there is an important 
exception to this result concerning the grid codes: the density- 
weighted velocity power spectrum produced by the grid code 
ENZO is almost identical to the power spectra produced with 
the SPH codes at f = 0.06 f cross , while FLASH, TVD, and ZEUS 
have lost a considerable amount of their power at k > 20. 
The power spectrum obtained with the ZEUS code shows the 
break into the dissipation range already at k ~ 10 and pro- 
duces a slightly steeper slope of about 1.65 than all other 
codes. At later times (f = 0.3 1 fcross-. and f = 0.62/ cross ) 
all codes produced similar density-weighted power spectra for 
k < 20, while the ENZO code develops a clear bottleneck (see e.g. 
Dobl er et all2003t lHaugen & Brandenburgl2004t ISchmidt et all 
2006a), which manifests itself in the excess power seen at k > 
10. Since E NZO was run here with a PPM diffusion parameter set 
to K = 0.0 dColella & Woodwardll 19841) the bottleneck effect i s 
quite strong (see also lKritsuk et al.l2007l: IFederrafh et alj20 09a). 
Although the FLASH code uses the same numerical technique as 
ENZO it does not show such a pronounced bottleneck effect, be- 
cause FLASH was use d with the recommended PPM diffusion 
parameter of K = 0.1 dColella & Woodward 19841) . 

As the decay progresses (f = 3.1, and 6.2 fcross); power gets 
dissipated differently by the various codes at k > 8. However, 
on large scales (k < 8), all codes used in the present study show 
very similar volume- and density-weighted velocity power spec- 
tra with slight variations that can be attributed to sta tistical fluc- 
tuations dKritsuk et alj|2007t IFederrafh etai1l2009al) . This is an 
important result, as it shows that all codes, despite having differ- 
ent dissipation mechanisms acting on small scales at k > 8, on 
large spatial scales the results of the decaying turbulence experi- 
ments presented here are quite robust. They do not show consid- 
erable systematic differences for the codes employed here at the 
resolutions studied. Moreover, it is important to note that the dis- 
sipation ranges at k > 8 are not just different when we compare 
grid codes with SPH codes, but they are also different among the 
SPH codes and among the grid codes. Thus, we conclude that the 
dissipation range is strongly dependent on the code being used, 
while scales with k < 8 are similarly well reproduced by all the 
hydrodynamic codes employed here. 

The overall performance of the codes, as seen through the 
analysis of their velocity power spectra at times f = 3.1f C ross, 
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and t = 6.2 f cross shows that the numerical viscosity of grid 
codes is generally lower than that of SPH codes, and details of 
the method used will determine the detailed ranking. For exam- 
ple, using K — in PPM, as done for the ENZO run, yields a 
lower dissipation value, although with a stronger bottleneck ef- 
fect than the K = 1 value used for the FLA SH run. Similarly, we 
find th at the viscosity implementation by iMorris & Monaghanl 
d 1997b used in the PHANTOM B and VINE runs is superior to that 
of GADGET and the PHANTOM A run. 

Considering the resolution of A^ 3 = 256 3 cells used in the 
present study (215 3 SPH particles interpolated to a 256 3 grid), 
the fact that the velocity power spectra are different at k > 8 
implies that one should be cautious with the interpretation of 
results obtained with power spectra at wavenumbers k > N/32. 
This means that length scales smaller than about 32 grid cells for 
grid codes (and SPH codes using a similar number of resolution 
elements interpolated to a grid of equivalent size), are affected by 
the individual dissipation mechanisms acting in hydrodynamical 
codes. In contrast, the results of the various codes are robust for 
k < N/32. This is encouraging, because the results of all the hy- 
drodynamical codes used here agree well in this regime, and one 
is free to choose a code for modelling supersonic turbulence as 
long as only results for scales k<N /32 are considered. However, 
this also means that one needs resolutions of at least 1024 3 grid 
cells to obtain roughly one full decade in length scale over which 
a power law could be fitted to turbulent velocity spectra. In prac- 
tice, this range turns out to be even smaller than one decade in 
length scales at a resolution of 1024 3 grid ce lls dKritsuk et alj 
2007; iKlein et al]l2007t iFederrath et aT]5009ah . 



5.2. Kinetic energy dissipation rates 

Tables [TJ and [2] list the integrals over the volume- and density- 
weighted velocity spectra, respectively for each code and run 
at each of the snapshots presented here. Up to t — 0.31 f C ross> 
all SPH codes dissipate volume-weighted power slightly faster 
than the grid codes, while for the density-weighted power, both 
SPH and grid codes dissipate kinetic energy at roughly the same 
rate. At t — 0.62 f croS s, however, ZEUS has dissipated about 15% 
more power in velocity fluctuations than the other codes, while 
all other grid codes have still dissipated less than the SPH codes. 
The density-weighted integral gives more similar results for all 
codes at all times analysed. At t = 3.1, and t = 6.2f cross , all 
codes have roughly dissipated the same amount of volume- and 
density-weighted power, except for the PHANTOM B run having 
kinetic energies about 16% larger than all other codes. 

5.3. Time evolution of the RMS Mach number 

The time evolution of the RMS Mach number is shown in 
Figure [9] The dotted line shows the expected power-law de- 
cay rate M rms oc t~y 2 for supersonic turbu lence (Mac Low et all 
1998; IStone et all [19981 iMac Lowl [19991) . starting at an RMS 
Mach number of ~ 3.9: 



M rms (t) = 3.9 | — + 1 
i 1, 



-1/2 



(14) 



Clearly, the RMS turbulent flow remains supersonic (i.e. Mrms > 
1) for all times analysed in the present study. However, turbu- 
lent velocity fluctuations become smaller on smaller scales. The 
transition scale separating supersonic motions on large scales 
and subs onic motions on small scales i s called the sonic scal e 
k s (e.g.. lVazquez-Semadeni et al .1 120031: IFederrath et alll2009ah . 



We computed an estimate of the sonic scale using the definition 
of equation (0. 

5.4. Time evolution of the sonic scale 

Table [3] lists the evolution of the sonic scale for all codes/runs 
and all snapshots. The sonic scale decreases fastest for 
PHANTOM A, PHANTOM B and ZEUS, and slowest for ENZO. 
GADGET, VINE, FLASH and TVD show quite similar results for 
the evolution of the sonic scale. The sonic scale does not differ 
considerably among the various codes at later times (f = 3. 1 f cross 
and t = 6.2 fcross). However, during the initial stages of the de- 
cay, there are differences up to 30%. This is partly a result of 
our computation of the sonic scale (cf. § 13.21 ). Since the differ- 
ent codes have quite different dissipation properties, the sonic 
scale is affected accordingly (cf. § 15. 11 1. Also note that the sonic 
scale is given as an integer. Thus, the fact that, for instance, ENZO 
maintains k s = 10 until t = 0.62 f closs does not imply that its sonic 
scale stays exactly the same for all times t < 0.62 f cross , but will 
have also decreased slightly. However, our grid of wavenumbers 
is binned such that only integer values of k are permitted, and 
thus, rounding errors introduce uncertainties of about 10% in 
the £ s -values at early times of the decay (f < lf cr0 ss)- 

5.5. Probability distribution functions 

5.5.1 . Probability distribution functions of the gas density 

Figure[l0]shows volume-weighted probability distribution func- 
tions (PDFs) of the gas density. Each panel shows the compar- 
ison of the density PDFs for all codes at t = 0.0, 0.06, 0.31, 
0.62, 3.1, and 6.2 f cross after they have been interpolated to grids 
of 256 3 cells. The density PDFs were computed from the loga- 
rithm of the density s = In (p/po)- The PDF p(s) is expected to 
follow roughly a Gaussian distribution 



p(s) ds - 



1 



exp 



(s - so) 2 



2<rl 



ds 



(15) 



mear 
1997 


value of s (e.g., [Vazquez-Semadeni 1994; Padoan et al. 
; IStone et al.1 11998: Mac Lowl 19991 iNordlund & Padoad 


1999 


: Ostriker et al. 19991 iKlessenl 


2000; Ostriker et al. 2001; 



Glover & Mac Lowil2007tlKritsuk et al . 2007; B eetz et alJl2008t 
Federrath et al.l l2008bt iLemaster & Stone) 120081: ISchmidt et al.1 
2009; Federrath et alj|2009at iGlover et all2009l 

The density PDFs of all codes show little variation around 
the peak of the distribution and at the high-density tail, and they 
are all roughly consistent with log-normal distributions, equa- 
tion (JT3J. The low density tails show stronger variations. This is 
because the low density tail is subject to str onger temporal vari- 
ations caused by inter mittent fluctuations (iRritsuk et al.l [2007; 
IFederrath e t al. 2009a). In Table we list the values of sq and 
(Tj for each code. Note that for a log-normal volume-weighted 
density distribution, so = _ o"?/2. The means and standard de- 
viations of the PDFs are similar for all codes and vary only by 
about 10%, except for our ZEUS run at t — 0.62 f cross , which has 
a mean value |sol about 28% larger than the average over all runs 
at that time. This appears slightly too high a variation to be at- 
tributed to temporal fluctuations. The difference of the density 
PDF around its peak obtained with the ZEUS run at t — 0.62 ? cross 
is also visible in Figure [10] (bottom left panel). However, at later 
times, the ZEUS density PDFs are almost identical to the ones 
obtained with the other codes. 
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One would expect that SPH codes can resolve high-density 
regions better than grid codes. This is because the effective SPH 
resolution increases with increasing density (cf. § 14.2.3b . On the 
other hand, low-density gas will become less resolved as parti- 
cles move towards higher density regions and away from low- 
density regions. However, the density PDFs of grid codes and 
SPH codes agree very well after interpolation of the particle 
data to grids. This may be caused by our interpolation proce- 
dure. We therefore would like to test the density PDF obtained 
from the SPH particle density directly without interpolation to 
a grid. The density PDF of SPH particles is naturally mass- 
weighted. Therefore, in order to obtain the volume-weighted 
density PDF directly from the particles we must weight the con- 
tribution of each particle into a density bin by the inverse of its 
density (for equal mass particles). Note that in general, volume- 
weighted PDFs, p v , are related to mass-weighted PD Fs, p m , by 
Pv = Pm/P (e.g. bstriker et al.ll200lt iLi et al.ll2003l) . This con- 
version is necessary as particles and grid cells do not sample the 
computational volume in the same manner, i.e. the particles are 
concentrated in high-density regions, whereas the grid cells sam- 
ple the volume homogeneously. The same applies for the distri- 
bution of Lagrangian tracer particles used in g rid simulations 
dFederrath et al.ll2008at IPrice & Federra th 2009). Therefore, we 
must compensate the particle PDF for the fact that most particles 
are located in high-density regions. 

Figure [TT] shows the density PDF of the GADGET run as ob- 
tained both after interpolation to a grid (black line), and directly 
obtained from the SPH particle density (red line). Additionally, 
we show the SPH density PDF obtained from the VINE run for 
comparison. The SPH density PDFs for GADGET and VINE were 
transformed into volume-weighted PDFs using p v = p,„/p to 
allow for a direct comparison of the grid and the SPH density 
PDFs. The grid-interpolated and the SPH density PDFs show no 
significant differences: the grid-interpolated PDF exhibits less 
scatter in the low-density regime than the SPH particle PDF, in- 
dicating that the grid-interpolated PDF samples the low-density 
regime slightly better than the particle-based PDF. On the other 
hand, the particle-based density PDF shows better sampling of 
the high-density tail and extends to slightly higher densities and 
lower probability densities. This is to be expected because the 
effective resolution is larger for the particle-based distribution at 
high densities, and SPH density peaks are smoothed to slightly 
smaller densities by the interpolation technique, equation (Q~|i. 
This analysis offers an additional illustration of the fact that the 
adaptivity of the SPH method is suppressed when the SPH infor- 
mation is interpolated onto a grid of resolution smaller than the 
effective SPH resolution (cf. § |4~2~3l 

It is an encouraging result that all codes are able to reproduce 
the general form of the density PDF quite well. This is important, 
because the turbulent density PDF is an essential ingredient to 
many models of star formation: to understa nd the mass distribu- 
tion of cores in molecular clouds and stars dPadoan & Nordlundl 
l2002[:lHennebelle&Chabrierll2008ll2009l). the star form ation 
rate dKrumholz & McKedl2005l: iPadoan & Nordlundll2009l) . the 
star formation efficiency (Elmegreen 2008), and the Kennicutt- 
Schmidt relation o n galactic scales ( Elmegreenll2002t iKravtsovl 
120031: lTassisll2007l) . 



5.5.2. Probability distribution functions of the rate of strain 

Figure [12]presents the trace free rate of strain PDFs computed as 
explained in § 13.41 Each panel shows the comparison of the PDFs 
of all codes atf = 0.0, 0.06, 0.31, 0.62, 3.1, and6.2 f cross . Within 
the first f cross (top row and bottom-left panel of Fig. [T2l . all SPH 



codes and ZEUS have PDFs that are slightly narrower than those 
of the remaining grid codes (top-right and bottom-left panels of 
Fig. [T2l . At later times (bottom-middle and bottom-right pan- 
els of Fig. [T2l . all codes appear to have distributions with sim- 
ilar widths, but with the PDFs of the SPH codes and the ZEUS 
code peaking at almost half the \S* \ value of the other grid codes. 
In particular, GADGET, PHANTOM A, PHANTOM B, VINE and ZEUS 
have their peaks at log|5""| < 0, while ENZO and FLASH have 
their peaks at log |5"*| > 0. TVD appears to peak in between. At 
t = 6.2 f cross the velocity power spectra of the codes are main- 
tained up to fc-values that can be ordered as follows (from higher 
to lower values): ENZO, FLASH, TVD, ZEUS, VINE, PHANTOM B, 
GADGET, PHANTOM A (cf. the bottom-right panels of Fig. and [8] 
in § 15. 11 1. This is the order of the |5*|-value of the peaks of the 
rate of strain PDFs at this time (bottom-right panel of Fig.[T2l. 

5.5.3. Probability distribution functions of the vorticity 

In Figure [13] we present the vorticity PDFs. The vorticity was 
computed as explained in § 13.41 Each panel shows the compari- 
son of the vorticity PDFs of all codes at t = 0.0, 0.06, 0.3 1 , 0.62, 
3.1, and 6.2 f closs . The vorticity and the trace free rate of strain 
PDFs show a similar behaviour: within the first f cross (top row 
and bottom-left panel of Fig.fTSl. the SPH codes and ZEUS show 
narrower distributions than the remaining grid codes. At later 
times (bottom-middle and bottom-right panels of Fig. [T3l ), the 
vorticity distributions of the SPH codes and ZEUS peak at almost 
half the vorticity value of the other grid codes. Again, the peaks 
of PHANTOM B, VINE, TVD, and ZEUS are bracketed by the peaks 
of the remaining codes, with GADGET and PHANTOM A peaking at 
the lowest values of the vorticity, and ENZO and FLASH peaking 
at the highest vorticity. 

Since the vorticity is related to the ability of the codes to 
model turbulent eddies, these comparisons indicate that SPH 
codes exhaust their ability earlier than grid codes, most likely be- 
cause of their excess viscosity acting on the smallest of these ed- 
dies and erasing them. As in the case of the rate of strain above, 
the order of the vorticity values of the peaks of the vorticity PDFs 
is the same as the order with which the codes maintain their ve- 
locity power spectra in the high-A: regime (cf. the bottom right 
panel of Fig . [7] and [8] and the bottom-right panel of Fig . [T"2"b . 

5.5.4. Probability distribution functions of the velocity 

Figure [T4l shows the PDFs of the velocity component v z . Each 
panel shows the comparison of the PDFs of all codes at times 
t = 0.0, 0.06, 0.31, 0.62, 3.1, 6.2f cross . For all these snapshots, 
the velocity PDFs are very similar for all codes. Apart from sta- 
tistical fluctuations showing up in the wings of the distributions, 
the velocity PDFs are roughly Gaussian distributions. As the 
standard deviation of the density PDFs (cf. Fig. ITOb , the stan- 
dard deviation of the velocity PDFs decreases with time, which 
simply reflects the turbulence decay (cf. Fig.[9]l. 



6. Computational efficiency of the codes 

Table [5] provides a summary of the computational efficiency of 
our codes/runs for the present setup. We remind the reader that in 
the current study all SPH codes used 215 3 resolution elements, 
while the grid codes used 256 3 resolution elements. We com- 
pensated for this discrepancy in the number of resolution ele- 
ments, as well as compensated roughly for the different CPU 
clock rates in the last row of Table [5] However, the various runs 
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were performed on different parallel machines throughout the 
world, some with, others without optimisations. Furthermore, 
different supercomputers use different hardware solutions for the 
parallelisation. Thus, the given numbers should only be taken 
as a rough estimate that may be accurate to within factors of a 
few. However, we emphasise that we have performed both the 
GADGET and the ENZO run on exactly the same supercomputing 
platform with the same optimisations, such that we can compare 
the performance of these two runs directly. 

The fastest of all runs was performed with the TVD grid code. 
It is roughly ten times faster than our ENZO run, about 15 times 
faster than FLASH, and about 27 times faster than the version of 
ZEUS used here. However, the official ZEUS-MP version is ex- 
pected to be faster and to scale better. Due to its specific de- 
sign and implementation, PHANTOM is the fastest of all the SPH 
codes employed, but still roughly 50% slower than the slowest 
grid code (our version of ZEUS). GADGET and VINE are about 16 
times, and about 30 times, respectively slower than PHANTOM. 
However, it must be remembered that the extra cost in SPH re- 
flects the fact that resolution elements are placed to follow the 
mass, and thus preferentially to resolve high-density regions. 
Thus, additional information is calculated on small scales. This 
information however does not enter this comparison as the parti- 
cles are interpolated onto a fixed grid. The extra cost for the SPH 
runs is similar to what ca n occur with grid-based Adaptive Mesh 
Refinement (AMR) (e.g. Ber ger & ColellallT9 89). because of the 
additional overhead to store a nd iterate over the AMR hierarchy 
(see e.g. ISchmidt et al.ll200 9). However, a quantitative analysis 
of the performance of AMR versus SPH is beyond the scope of 
this paper, and should be discussed elsewhere. 

The fact that GADGET and VINE are more than one order 
of magnitude slower than any of the grid codes makes it com- 
putationally expensive to study supersonic turbulence at reso- 
lutions higher than about 256 3 SPH particles. This may partly 
explain why no SPH calculations of supersonic turbulence have 
so far been using more than 512 3 part icles. The latter w a s only 
achieved with the PHANTOM code in iPrice & Federrathl d2009l) 
and in the KITP code comparison projecfj 

7. Conclusions 

In this paper we report the comparison of the performance of 
four grid-based and three particle-based hydrodynamic codes on 
the modelling of supersonic turbulence decay. In particular, we 
have studied the decay of compressible, supersonic, isothermal 
turbulence in the absence of gravity within a periodic box using 
simulations with resolutions of 256 3 grid cells for the grid codes, 
and 215 3 particles for the SPH codes. 

We have simulated driven turbulence with the SPH code 
GADGET. The SPH particle distribution at the end of the driving 
has been interpolated onto a grid that provides the initial con- 
ditions to the grid codes employed, namely ENZO, FLASH, TVD 
and ZEUS. We have also followed the decay of turbulence using 
several implementations of the SPH method, namely GADGET, 
PHANTOM (runs: A, B), and VINE, where PHANTOM B used the 
M orris & Monaghanl d 19971) viscosity switch, while PHANTOM A 
was run without the switch. 

The turbulence decay was followed for about six turbulent 
crossing times. During the whole decay phase considered here, 
the turbulent flow stays in the supersonic regime. The turbulent 
energy dissipation was measured for all codes. For times greater 

6 http : //kitpstarformation07 . wiki spaces . com/ 
Star+Formation+Test+Problems 



than about 3 crossing times (when any initial transient phase due 
to small variations in the initial conditions has disappeared, and 
all codes have developed their individual dissipation signatures) 
a comparison of volume- and density-weighted velocity spectra 
indicates that the numerical viscosity of grid codes is generally 
lower than that of SPH codes, with details of the method like 
the order of the code contributing secondarily. We show that the 
differences between ENZO and FLASH are due to our choice of 
using different PPM diffusion parameters for the two codes, i.e., 
in this study ENZO was used with a PPM diffusion parameter of 
K = 0.0, while FLASH was used with the PP M diffusion param- 
eter se t to K — 0.1 as recommended by Colella & Woodward] 
(119841) . Switching-off PPM diffusion completely (as in our ENZO 
run) results in less dissipation, but produces a stronger bottle - 
neck effec t (see alsolKritsuk et alJl2007tlFederrath et aLl l2009a). 
Use of the Morris & Monaghan ( 1997) viscosity implementation 
for SPH provides less dissipation as observed in our PHANTOM B 
and VINE runs in comparison with the GADGET and PHANTOM A 
runs. In general, the viscosity acts differently for the different 
grid- and SPH-codes at wavenumbers k > 8 at the resolutions 
studied here, shown by our analysis of velocity power spectra in 
§ 15.11 However, all codes produced velocity spectra that are in 
good agreement for k < 8. 

Using Fourier spectra we also showed that the additional in- 
formation that the SPH method can offer in high-density regions 
and/or on small scales will be suppressed if it is not interpolated 
onto a high enough resolution grid, as discussed in § 14.2.31 and 

§o 

The trace-free rate of strain and the vorticity PDFs confirm 
the ordering of the runs according to their dissipation given 
above. The density PDFs are very similar for all the runs per- 
formed in the present study. The means and standard deviations 
of the logarithmic density varied by less than 10% for all codes 
at all times analysed, with one exception (the mean logarith- 
mic density obtained in our ZEUS run at t — 0.62 f cross varied 
by about 30%, which is also seen in its density PDF). For the 
SPH code GADGET we have shown that the density PDF obtained 
from the SPH particle distribution samples the high density tail 
slightly better, while our results indicate that using a grid, the 
low-density tail is slightly better sampled. However, the overall 
shape is very close to a log-normal distribution in density, and 
its mean and standard deviation are quite robust for all codes 
employed in the present study. 

Our results demonstrate that different codes have differ- 
ent dissipation mechanisms affecting spatial scales k > N/32. 
However, our code comparison also shows that SPH and grid 
codes give similar results for an equivalent number of resolu- 
tion elements N for each direction in space on scales k < N/32, 
though with the SPH runs being about ten times more computa- 
tionally expensive than the grid runs on average. Careful choice 
of numerical algorithm can extend this scaling range slightly, 
indicating that grid codes tend to show a slightly longer scal- 
ing range than SPH codes (cf. Figures|7]and[8]l. However, at the 
numerical resolutions employed in the present study, all slopes 
inferred from the volume- and density-weighted spectra are too 
steep compared with higher resolution data from the literature. 
It is thus rather a question of resolution than a question of the 
specific properties of the hydrodynamical codes used in this 
study that determines their performance in reproducing turbu- 
lence scaling relations. This must be tested in a future compari- 
son employing at least one order of magnitude more resolution 
elements for both grid and SPH codes. 

Acknowledgements. The authors would like to thank Jens Niemeyer for many 
interesting and stimulating discussions, especially during the setup phase of 



S. Kitsionas et al.: Algorithmic comparisons of decaying, isothermal, supersonic turbulence 



11 



this project. We thank the anonymous referee for a balanced and detailed re- 
port. SK kindly acknowledges support by an EU Commission "Marie Curie 
Intra-European (Individual) Fellowship" of the 6th Framework Programme with 
contract number MEIF-CT-2004-01 1226. SK also acknowledges financial as- 
sistance by the EU Commission Research Training Network "Constellation" 
of the 6th Framework Programme. CF acknowledges financial support from 
the International Max Planck Research School for Astronomy and Cosmic 
Physics (IMPRS-A) and the Heidelberg Graduate School of Fundamental 
Physics (HGSFP). The HGSFP is funded by the Excellence Initiative of the 
German Research Foundation DFG GSC 129/1. The ENZO and GADGET sim- 
ulations used computational resources from the HLRBII project grant h0972 
and pr321o at the Leibniz Rechenzentrum Garching and from a project grand 
by CASPUR, Italy. CF and RSK are grateful for subsidies from the DFG 
SFB 439 Galaxies in the Early Universe. RSK and CF acknowkledge finan- 
cial support from the German Bundesministerium fur Bildung und Forschung 
via the ASTRONET project STAR FORMAT (grant 05A09VHA) and from 
the Deutsche Forschungsgemeinschaft (DFG) under grants no. KL 1358/1, KL 
1358/4, KL 1359/5. RSK furthermore thank for subsidies from a Frontier grant 
of Heidelberg University sponsored by the German Excellence Initiative and 
for support from the Landesstiftung Baden- Wiirttemberg via their programme 
International Collaboration II. DJP is supported by a Royal Society University 
Research Fe llowship (UK). The PDF figures were partly produced using splash 
iPrice 2007). MG and SW acknowledge support by the DFG cluster of ex- 
cellence "Origin and Structure of the Universe". All VINE and part of the 
FLASH calculations were performed on an SGI Altix 3700 Bx2 supercom- 
puter that was funded by the DFG cluster of excellence "Origin and Structure 
of the Universe". The simulations performed by JK utilized a high perfor- 
mance cluster that was built with funding from the Korea Astronomy and Space 
Science Institute (KASI) and the Korea Science and Engineering Foundation 
through the Astrophysical Research Center for the Structure and Evolution of 
Cosmos (ARCSEC). The work of JK was supported by ARCSEC. AKJ ac- 
knowledges support by the Human Resources and Mobility Programme of the 
European Community under the contract MEIF-CT-2006-039569. MMML ac- 
knowledges partial support for his work from NASA Origins of Solar Systems 
grant NNX07AI74G. The software used in this work was in part developed by 
the DOE-supported ASC / Alliance Center for Astrophysical Thermonuclear 
Flashes at the University of Chicago. 



References 

Attwood, R. E., Goodwin, S. P., & Whitworth, A. P. 2007, A&A, 464, 447 
Ballesteros-Paredes, J., Gazol, A., Kim, J., et al. 2006, ApJ, 637, 384 
Ballesteros-Paredes, J., Klessen, R. S., Mac Low, M.-M., & Vazquez-Semadeni, 

E. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 

63-80 

Beetz, C, Schwarz, C, Dreher, J., & Grauer, R. 2008, Physics Letters A, 372, 
3037 

Benzi, R., Biferale, L., Fisher, R. T, et al. 2008, Physical Review Letters, 100, 
234503 

Berger, M. J. & Colella, P. 1989, Journal of Computational Physics, 82, 64 
Blitz, L., Fukui, Y., Kawamura, A., et al. 2007, in Protostars and Planets V, ed. 

B. Reipurth, D. Jewitt, & K. Keil, 81-96 
Boldyrev, S. 2002, ApJ, 569, 841 

Boldyrev, S., Nordlund, A., & Padoan, P. 2002a, ApJ, 573, 678 
Boldyrev, S., Nordlund, A., & Padoan, P. 2002b, Physical Review Letters, 89, 
031102 

Bonazzola, S., Perault, M., Puget, J. L., et al. 1992, Journal of Fluid Mechanics, 
245, 1 

Calder, A. C, Fryxell, B., Plewa, T., et al. 2002, ApJS, 143, 201 
Chandrasekhar, S. 1949, ApJ, 110, 329 

Chandrasekhar, S. 1951a, Royal Society of London Proceedings Series A, 210, 
18 

Chandrasekhar, S. 1951b, Royal Society of London Proceedings Series A, 210, 
26 

Colella, P. & Woodward, P. R. 1984, Journal of Computational Physics, 54, 174 
Commercon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2008, 

A&A, 482, 371 
de Avillez, M. A. & Mac Low, M.-M. 2002, ApJ, 581, 1047 
Dimonte, G., Youngs, D. L., Dimits, A., et al. 2004, Physics of Fluids, 16, 1668 
Dobler, W., Haugen, N. E. L., Yousef, T. A., & Brandenburg, A. 2003, Phys. Rev. 

E, 68, 026304 

Dubey, A., Fisher, R., Graziani, C, et al. 2008, in Astronomical Society of the 
Pacific Conference Series, Vol. 385, Numerical Modeling of Space Plasma 
Flows, ed. N. V. Pogorelov, E. Audit, & G. P. Zank, 145 

Elmegreen, B. G. 2002, ApJ, 577, 206 

Elmegreen, B. G. 2008, ApJ, 672, 1006 

Elmegreen, B. G. & Scalo, J. 2004, ARA&A, 42, 21 1 



Federrath, C, Duval, J., Klessen, R. S., Schmidt, W., & Mac Low, M.-M. 2009a, 

A&A, submitted (arXiv: 0905. 1060 1 
Federrath, C, Glover, S. C. O., Klessen, R. S., & Schmidt, W. 2008a, Phys. Scr. 

T, 132,014025 

Federrath, C, Klessen, R. S., & Schmidt, W. 2008b, ApJ, 688, L79 

Federrath, C, Klessen, R. S., & Schmidt, W. 2009b, ApJ, 692, 364 

Frisch, U. 1995, Turbulence (Cambridge Univ. Press) 

Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273 

Gingold, R. A. & Monaghan, J. J. 1977, MNRAS, 181, 375 

Glover, S. C. O., Feder rath, C, Mac L ow, M.-M., & Klessen, R. S. 2009, 

MNRAS, accepted larXiv:0907.4081 1 
Glover, S. C. O. & Mac Low, M.-M. 2007, ApJ, 659, 1317 
Gritschneder, M., Naab, T., Walch, S., Burkert, A., & Heitsch, F. 2009, ApJ, 694, 

L26 

Haugen, N. E. & Brandenburg, A. 2004, Phys. Rev. E, 70, 026405 
Hayes, J. C, Norman, M. L., Fiedler, R. A., et al. 2006, ApJS, 165, 188 
Heitmann, K, Ricker, P. M., Warren, M. S., & Habib, S. 2005, ApJS, 160, 28 
Heitsch, E, Mac Low, M.-M., & Klessen, R. S. 2001, ApJ, 547, 280 
Hennebelle, P. & Chabrier, G. 2008, ApJ, 684, 395 
Hennebelle, P. & Chabrier, G. 2009, ApJ, 702, 1428 
Isichenko, M. B. 1992, Reviews of Modern Physics, 64, 961 
Jappsen, A.-K., Klessen, R. S., Larson, R. B., Li, Y, & Mac Low, M.-M. 2005, 
A&A, 435, 611 

Kim, J., Ryu, D., Jones, T. W., & Hong, S. S. 1999, ApJ, 514, 506 

Kitsionas, S. & Whitworth, A. P. 2002, MNRAS, 330, 129 

Klein, R. I., Inutsuka, S.-L, Padoan, P., & Tomisaka, K. 2007, in Protostars and 

Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 99-1 16 
Klessen, R. S. 2000, ApJ, 535, 869 

Klessen, R. S., Heitsch, E, & Mac Low, M.-M. 2000, ApJ, 535, 887 
Klessen, R. S. & Lin, D. N. 2003, Phys. Rev. E, 67, 04631 1 
Kolmogorov, A. N. 1941, Dokl. Akad. Nauk SSSR, 32, 16 
Kowal, G. & Lazarian, A. 2007, ApJ, 666, L69 
Kravtsov, A. V. 2003, ApJ, 590, LI 

Kritsuk, A. G., Norman, M. L., & Padoan, P. 2006, ApJ, 638, L25 

Kritsuk, A. G., Norman, M. L., Padoan, P., & Wagner, R. 2007, ApJ, 665, 416 

Krumholz, M. R. & McKee, C. F. 2005, ApJ, 630, 250 

Lemaster, M. N. & Stone, J. M. 2008, ApJ, 682, L97 

Lemaster, M. N. & Stone, J. M. 2009, ApJ, 691, 1092 

Lesieur, M. 1997, Turbulence in Fluids (Kluwer) 

Li, Y, Klessen, R. S., & Mac Low, M.-M. 2003, ApJ, 592, 975 

Lucy, L. B. 1977, AJ, 82, 1013 

Mac Low, M.-M. 1999, ApJ, 524, 169 

Mac Low, M.-M. & Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125 
Mac Low, M.-M., Klessen, R. S., Burkert, A., & Smith, M. D. 1998, Physical 

Review Letters, 80, 2754 
McKee, C. F. & Ostriker, E. C. 2007, ARA&A, 45, 565 
Metzler, R. & Klafter, J. 2000, Phys. Rep., 339, 1 
Monaghan, J. J. 1988, Computer Physics Communications, 48, 89 
Monaghan, J. J. 1992, ARA&A, 30, 543 

Monaghan, J. J. 1997, Journal of Computational Physics, 136, 298 
Monaghan, J. J. 2005, Reports on Progress in Physics, 68, 1703 
Monaghan, J. J. & Gingold, R. A. 1983, Journal of Computational Physics, 52, 
374 

Monaghan, J. J. & Lattanzio, J. C. 1985, A&A, 149, 135 

Morris, J. & Monaghan, J. 1997, Journal of Computational Physics, 136, 41 

Nelson, A. F, Wetzstein, M., & Naab, T. 2008, arXiv:0802.4253 

Nordlund, A. & Padoan, P. 1999, in Interstellar Turbulence, ed. J. Franco & 

A. Carraminana, 218 
Norman, M. L. & Bryan, G. L. 1999, in Astrophysics and Space Science 

Library, Vol. 240, Numerical Astrophysics, ed. S. M. Miyama, K. Tomisaka, 

& T. Hanawa, 19 
O'Shea, B. W., Bryan, G., Bordner, J., et al. 2004, arXiv:0403044 
Ostriker, E. C, Gammie, C. E, & Stone, J. M. 1999, ApJ, 513, 259 
Ostriker, E. C, Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980 
Padoan, P., Jimenez, R., Nordlund, A., & Boldyrev, S. 2004, Physical Review 

Letters, 92, 191102 
Padoan, P. & Nordlund, A. 2002, ApJ, 576, 870 
Padoan, P. & Nordlund, A. 2009, ApJ, submitted larXiv:0907.0248 I 
Padoan, P., Nordlund, A., & Jones, B. J. T. 1997, MNRAS, 288, 145 
Padoan, P., Nordlund, A., Kritsuk, A. G., Norman, M. L., & Li, P. S. 2007, ApJ, 

661, 972 

Panis, J.-F. & Perault, M. 1998, Physics of Fluids, 10, 3111 
Porter, D. H. & Woodward, P. R. 1994, ApJS, 93, 309 
Price, D. J. 2004, PhD Thesis, University of Cambridge 
Price, D. J. 2007, Publ. Astron. Soc. Aust., 24, 159 
Price, D. J. & Federrath, C. 2009, MNRAS, submitted 
Price, D. J. & Monaghan, J. J. 2004, MNRAS, 348, 139 
Price, D. J. & Monaghan, J. J. 2007, MNRAS, 374, 1347 



12 



S. Kitsionas et al.: Algorithmic comparisons of decaying, isothermal, supersonic turbulence 



Sasao, T. 1973, PASJ, 25, 1 

Scalo, J. & Elmegreen, B. G. 2004, ARA&A, 42, 275 

Schmidt, W., Federrath, C, Hupp, M., Kern, S., & Niemeyer, J. C. 2009, A&A, 
494, 127 

Schmidt, W., Federrath, C, & Klessen, R. 2008, Phys. Rev. Lett., 101, 194505 
Schmidt, W., Hillebrandt, W., & Niemeyer, J. C. 2006a, Computers and Fluids, 
35, 353 

Schmidt, W., Niemeyer, J. C, & Hillebrandt, W. 2006b, A&A, 450, 265 
Schmidt, W., Niemeyer, J. C, Hillebrandt, W., & Ropke, F. K. 2006c, A&A, 450, 
283 

She, Z.-S. & Leveque, E. 1994, Physical Review Letters, 72, 336 

Springel, V., Yoshida, N., & White, S. D. M. 2001, New Astronomy, 6, 79 

Stone, J. M. & Norman, M. L. 1992a, ApJS, 80, 753 

Stone, J. M. & Norman, M. L. 1992b, ApJS, 80, 791 

Stone, J. M., Ostriker, E. C, & Gammie, C. F. 1998, ApJ, 508, L99 

Strang, G. 1968, SIAM J. Numer. Anal., 5, 506 

Tasker, E. J., Brunino, R., Mitchell, N. L., et al. 2008, MNRAS, 390, 1267 
Tassis, K. 2007, MNRAS, 382, 1317 

Toro, E. F. 1997, Riemann solvers and numerical methods for fluid dynamics 
(Springer) 

van Leer, B. 1977, Journal of Computational Physics, 23, 276 
Vazquez-Semadeni, E. 1994, ApJ, 423, 681 

Vazquez-Semadeni, E., Ballesteros-Paredes, J., & Klessen, R. S. 2003, ApJ, 585, 
L131 

Vernaleo, J. C. & Reynolds, C. S. 2006, ApJ, 645, 83 
von Weizsacker, C. F. 1951, ApJ, 114, 165 

von Weizsacker, C. F. V. 1943, Zeitschrift fur Astrophysik, 22, 319 
Wetzstein, M., Nelson, A. F, Naab, T., & Burkert, A. 2008, arXiv:0802.4245 



S. Kitsionas et al.: Algorithmic comparisons of decaying, isothermal, supersonic turbulence 
Table 1. Time evolution of E tot = J E(k) die in units of [c^], for all codes/runs. 
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Time t [f cross ] GADGET PHANTOM A PHANTOM B VINE ENZO FLASH TVD ZEUS 



0.0 7.64 7.64 7.64 7.64 7.64 7.64 7.64 7.64 

0.06 7.16 7.12 7.12 7.31 7.58 7.60 7.61 7.51 

0.31 5.72 5.72 5.71 5.83 6.32 6.16 6.17 5.88 

0.62 4.68 4.72 4.72 4.76 5.22 5.03 5.04 4.51 

3.1 1.72 1.70 1.74 1.71 1.75 1.74 1.75 1.64 

6.2 1.16 1.14 1.32 1.15 1.14 1.14 1.14 1.15 



Table 2. Time evolution of E mV[ ,^ ot = J E mw (k) dk in units of [po cj], for all codes/runs. 



Time t [t aoss ] GADGET PHANTOM A PHANTOM B VINE ENZO FLASH TVD ZEUS 



0.0 7.22 7.22 7.22 7.22 7.22 7.22 7.22 7.22 

0.06 6.84 6.93 6.95 7.00 7.00 6.82 6.83 6.82 

0.31 5.42 5.49 5.50 5.54 5.70 5.41 5.43 5.39 

0.62 4.30 4.36 4.38 4.38 4.55 4.31 4.32 4.17 

3.1 1.56 1.56 1.61 1.56 1.57 1.56 1.55 1.56 

6.2 1.12 1.11 1.29 1.12 1.11 1.10 1.10 1.13 



Table 3. Time evolution of the sonic scale k s as defined by equation (0 for all codes/runs. 



Time t [t aoS!} ] GADGET PHANTOM A PHANTOM B VINE ENZO FLASH TVD ZEUS 



0.0 10 10 10 10 10 10 10 10 

0.06 9 9 9 9 10 10 10 9 

0.31 9 8 8 8 10 998 

0.62 8 7 7 8 10 997 

3.1 3 3 3 33333 

6.2 2 2 2 22222 



Table 4. Time evolution of the mean sq and standard deviation cr s of the logarithmic density s — In (p/po) for all codes/runs. 



GADGET PHANTOM A PHANTOM B VINE ENZO FLASH TVD ZEUS 

Time t [f cross ] s <r s s a s s <r s s cr s s cr s s <r s s cr s s cr s 

0.0 -0.73 1.23 -0.73 1.23 -0.73 1.23 -0.73 1.23 -0.73 1.23 -0.73 1.23 -0.73 1.23 -0.73 1.23 

0.06 -0.72 1.22 -0.71 1.21 -0.71 1.21 -0.72 1.22 -0.73 1.23 -0.73 1.23 -0.73 1.23 -0.73 1.23 

0.31 -0.72 1.22 -0.70 1.20 -0.70 1.21 -0.72 1.22 -0.73 1.23 -0.73 1.23 -0.73 1.23 -0.73 1.23 

0.62 -0.54 1.10 -0.52 1.08 -0.53 1.08 -0.54 1.09 -0.58 1.14 -0.58 1.13 -0.59 1.14 -0.68 1.19 

3.1 -0.15 0.56 -0.15 0.55 -0.15 0.56 -0.15 0.56 -0.16 0.56 -0.16 0.56 -0.16 0.58 -0.16 0.58 

6.2 -0.06 0.34 -0.05 0.32 -0.06 0.33 -0.06 0.34 -0.06 0.34 -0.06 0.34 -0.06 0.35 -0.06 0.36 



Table 5. Computational efficiency of all codes/runs. 





GADGET 


PHANTOM A 


PHANTOM B 


VINE 


ENZO 


FLASH 


TVD 


ZEUS 


architecture (CPU type) 


Itanium 


Clovertown 


Clovertown 


Itanium 


Itanium 


Xeon 


Xeon 


Xeon 


CPU clock rate 


1.6 GHz 


2.66 GHz 


2.66 GHz 


1.3 GHz 


1.6 GHz 


1.6 GHz 


2.66 GHz 


2.33 GHz 


number of CPUs used 


32 


8 


8 


64 


8(32) 


64 


8 


64 


parallelisation 


MPI 


OpenMP 


OpenMP 


OpenMP 


MPI 


MPI 


MPI 


MPI 


total CPU-h 


6,490 


248 


248 


15,000 


165 (203) 


256 


10 


315 


total CPU-h (norm.") 


10,960 


697 


697 


20,600 


165 (203) 


256 


17 


459 



" CPU hours normalised to 256 3 resolution elements, and normalised to a clock rate of 1.6 GHz on an Intel Itanium/Xeon chip. These are very 
rough estimates that should only be accurate to within factors of a few, because of the different parallelisation hardware used on the various 
supercomputing platforms. Note however that GADGET and ENZO were run on the same supercomputing platform, and with the best run time and 
parameter optimisations that we could achieve on that supercomputer for both codes. The total CPU time for these two runs can thus be compared 
directly. ENZO was also run on 32 CPUs (values given in brackets). 
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Fig. 1. Time evolution of the RMS Mach number during the driving phase. After driving for about 2 turbulent crossing times 
fcross (see text), the RMS Mach number has reached a statistical steady state. The initial conditions for the decaying turbulence 
code comparison were chosen randomly at t — 2.5f cl0SS in the regime of fully developed supersonic turbulence, t > 2f cross (e.g. 
iFederrath et al.ll2009bl) . when the RMS Mach number has reached its statistical steady state of M lms ~ 3.9. 




Fig. 2. Top: Velocity power spectra obtained at four different snapshots during the driving phase (volume-weighted spectra on the 
left, and density-weighted spectra on the right panels). The spectrum of each snapshot is plotted with a different colour (1.2? cross : 
black lines; 2.5 f cross : red lines; 3.1 t aoss : green lines; 3.7 f cross : blue lines). The decay experiments using the SPH and grid codes 
starts with the snapshot at 2.5 ? cross when turbulence is fully established. The solid lines correspond to E(k) = E so i + £ com (volume- 
weighted, left), and E mw (k) = E mv/ , so i + £ m w,com (density-weighted, right), the dashed lines to the solenoidal (transverse) part, and 
the dash-dotted lines to the compressible (longitudinal) part of the spectra. Bottom: The same spectra as on the top, but compensated 
with power-law slopes of 2.20 (left panel - volume- weighted case) and 1.67 (right panel - density-weighted case). 




Fig. 3. Velocity power spectra of the initial conditions used for the decay experiments (volume- weighted spectrum on the left, 
and density-weighted spectrum on the right panel). The spectra were compensated with power-law slopes of 2.20 (left panel - 
volume-weighted case) and 1.67 (right panel - density- weighted case). 




Fig. 4. Left panel: Density-weighted velocity power spectra of the initial conditions using different velocity weights: v, 
(p/Po)'^ 2 v (black lines), and v mw = (p/Po)^ 3 v ( re d lines). Right panel: Density-weighted velocity power spectra [v, 
(p/po) 1 ^ 3 v] of the initial conditions interpolated to grids of 256 3 cells (black lines) and 512 3 cells (red lines). 
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Fig. 5. Time evolution of the volume-weighted velocity power spectra (compensated with power-law slopes of 2.20) of the SPH 
codes GADGET, PHANTOM (runs: A, B), and VINE (left column) and of the grid codes ENZO, FLASH, TVD, and ZEUS (right column) for 
the following snapshots: t = 0.0, 0.06, 0.31, 0.62, 3.1, and 6.2 f cr0 ss- 
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Fig. 6. Same as Fig. [5] but the density-weighted velocity power spectra (compensated with power-law slopes of 1.67) are shown. 
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Fig. 7. Comparison of the volume-weighted velocity power spectra (compensated with power-law slopes of 2.20) of all codes 
different times along the decay: t = 0.0, 0.06, 0.31, 0.62, 3.1, and 6.2 ? croS s- 
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Fig. 8. Same as Fig. [7] but the density- weighted velocity power spectra (compensated with power-law slopes of 1.67) are shown. 
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Fig. 9. Evolution of the RMS Mach number as a function of time in units of the turbulent cros sing time f cross for all codes/runs. 
The dotted line show s the expected power-law decay rate M rms oc r 1 ^ 2 for supersonic turbulence (Ma c Low et al.llT998b IStone et al.l 
1998: lMacLowlll999h . 
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Fig. 10. Comparison of the volume-weighted density PDFs of all codes at different times along the decay: t = 0.0, 0.06, 0.31, 0.62, 
3.1, and 6.2 ? cross . 
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Fig. 11. Volume-weighted density PDFs at t — 0.0, 0.06, 0.31, 0.62, 3.1, and 6.2 f cross based on all cells of the interpolated grid 
(black line), and based on all SPH particle densities for GADGET (red line), and all SPH particle densities for VINE (blue line). The 
SPH particle density PDFs were computed on the SPH particles directly. Thus, they do not involve any interpolation to a grid, but 
they needed to be converted to a volume-weighted density PDF by using p v = p m /p, in order to allow for a direct comparison to the 
grid-interpolated PDF (black line). Both the density PDF computed from the grid and the SPH densities agree very well, with the 
SPH density PDF extending to slightly higher densities, and exhibiting slightly less scatter at the high-density tail. 
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12. Same as Figure[lOl but the PDFs of the trace-free rate of strain are shown, as defined in eq. ( fT3l >. 
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13. Same as Figure[TfjJ but the PDFs of the vorticity are shown. 



